Quantum Monte Carlo method for models of molecular nanodevices 
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We introduce a quantum Monte Carlo technique to calculate exactly at finite temperatures the 
Green function of a fermionic quantum impurity coupled to a bosonic field. While the algorithm 
is general, we focus on the single impurity Anderson model coupled to a Holstein phonon as a 
schematic model for a molecular transistor. We compute the density of states at the impurity in a 
large range of parameters, to demonstrate the accuracy and efficiency of the method. We also obtain 
the conductance of the impurity model and analyze different regimes. The results show that even 
in the case when the effective attractive phonon interaction is larger than the Coulomb repulsion, a 
Kondo-like conductance behavior might be observed. 
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The physics of quantum impurities is attracting a great 
deal of attention. For decades, their study has been a 
classic subject in the area of strongly correlated electron 
systems Q. However, with the dramatic recent advances 
in nanoscale physics their interest has become paramount 
to a much broader audience. Among the most interest- 
ing recent achievements in nanometre-scale devices we 
find the ability to attach individual molecules to metallic 
electrodes 2J. In these systems, the coupling between 
the discrete internal degrees of freedom of the molecule 
and the macroscopic leads gives rise to new transport 
channels enabled by resonances originated in many-body 
interactions. An important example is the observation of 
the physics of the Kondo effect Q , which were previously 
seen in quantum dots |4j . In general, one expects that the 
transport through molecular devices will be affected by 
the coupling between their internal vibration modes and 
the electrons 0- In addition, one should also consider 
the possible electron correlation effects due to Coulomb 
repulsion within the molecule 

The solution of quantum impurity models has re- 
mained a difficult task and exact solutions are only 
possible in very specific instances with the use of tech- 
niques such as the Bethe ansatz. For general models, 
however, theoreticians have to relay heavily on numerical 
methods. Among those we count the celebrated Wilson's 
numerical renormalization group (NRG) and its various 
generalizations 0. This method is extremely accurate 
for the determination of the low energy physics of the 
model, but has limitations for the investigation of inter- 
mediate temperatures and impurity models with many 
orbitals. On the other hand, there is a very efficient 
Quantum Monte Carlo (QMC) method introduced by 
Hirsch and Fye 5\ which does not present those limita- 
tions and neither the so called "negative sign" problem. 

Among models of quantum impurities a most impor- 
tant one is the single impurity Anderson model (SIAM) 
Hamiltonian , that has been crucial for a wide number 



of physical problems. Essentially it consists of an atomic 
orbital (the impurity) with an on-site repulsive Coulomb 
interaction that is hybridized with a band of conduction 
electrons. 

The goal of the present work is twofold: we shall 
first introduce a new algorithm which generalizes that 
of Hirsch and Fye to the quantum impurity problem of 
a SIAM coupled to a general bosonic field. In order to 
test our method we shall consider the particular case of 
the SIAM coupled to a Holstein phonon which has been 
recently studied with NRG at T — 0. We shall present 
the detailed behavior of the density of states (DOS) of 
the impurity for various phonon strengths, local Coulomb 
interaction and temperatures. The second goal of our 
work will be to discuss new results for the conductance 
of this quantum impurity as a function of the gate volt- 
age, that can be understood under the light of our DOS 
calculations. We shall show the two generic conductance 
behaviors that occur depending on the values of the pa- 
rameters. Our results should be useful for interpretation 
of experiments of transport through molecular transis- 
tors. 

The SIAM coupled to a general bosonic mode is de- 
scribed by the action: 

S = I ' dTdr%{r)G^ l {T-r')c a {T') 

+ J dTdT<j)(T)Dv 1 (T-T')<j)(T')+ J dTUrif(T)ni(T) 

+ Xj2fdr4>(T)(c a (r)cAT)-\), (1) 

where c CT , c a are Grassmann variables for the creation and 
destruction of an electron with spin a at the impurity and 
cj> is the bosonic field. The coupling A is the strength of 
the interaction between the electrons at the impurity and 
the bosonic degree of freedom, and U is the Coulomb cost 
for the double occupation of the impurity site. The non- 
interacting propagators are Do(u> m ) for the bosons and 
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Go(u} n ) = [iiv n — V g + A(w„)] _1 for the impurity electrons, 
being u m and tv n bosonic and fermionic Matsubara fre- 
quencies respectively. The hybridization function A(uj n ) 
describes the coupling to the (given) conduction band 
and V g denotes the gate voltage. 

We shall now describe the QMC algorithm. As usual 
one performs a Trotter breakup of imaginary time into 
M time slices t; of length At — (3/M to discretize the ac- 
tion m and introduce a discrete Hubbard-Stratonovich 
transformation with auxiliary Ising-like fields s; , I = 
0, ...,M- 1. 

This renders the action quadratic in the fermionic vari- 
ables which can be integrated out. The resulting parti- 
tion function is 

Z = Tr {si} J P0Dei[Gf^^)]Dei[G^(T;^')] 

eXP {~R ^^D^iiuJm^m + A0 O }, (2) 

m 

where <fi m are the Fourier transform components of the 
bosonic field that obey 4> m — 4>_ m - The s; and <p m de- 
pendent electron Green functions are 

G *«A(n,n>) = Go X (r; - n>) + 5 u ,[aX uSl + AtA0(t ; )], 

(3) 

with coshAu = exp (AtU /2). A crucial new aspect of our 
algorithm is the choice of a simultaneous representation 
in terms of the Ising-like and continuous bosonic fields. 
Other algorithms chose to integrate out the bosonic vari- 
ables to obtain a time dependent electron density-density 
interaction which is then decoupled via a single contin- 
uous Hubbard-Stratonovich field. That procedure gives 
up the very nice properties of the discrete Ising-like fields 
and should be expected to have serious autocorrelation 
problems. In the present case, each degree of freedom is 
simulated in the most effective way. This is evident for 
the Ising-fields, but should also be clear for the bosonic 
fields as our frequency modes representation is reminis- 
cent of the Fourier acceleration method. 'In fact, while 
the update of the Ising fields is known to lead to very 
short autocorrelation times the same is not true for the 
bosonic fields. The physical reason is that the most rele- 
vant bosonic configurations are those in which </>(t) have 
a slow variation with r. To generate those configurations 
one needs to propose moves that coherently and unbiased 
involve all the M variables 4>{ti). Our algorithm provides 
a systematic way to do this by updating the frequency 
components 4> m of the bosonic fields. It is important 
to mention that the autocorrelation times of cf> m remain 
pretty long (typically about a hundred sweeps) but, un- 
like what would be updating in time domain, practically 
tractable.' 

The integrand of (J2J) thus defines the Boltzmann weight 
yV[si;4> m ] for the stochastic evaluation of the trace on 
the Ising and bosonic fields. We use a heat-bath algo- 
rithm according to which new configurations {si;cj) m } 



are generated and acepted with a probability P(W — > 
W) = W/(W + W). The ensuing Markov chain is: 

(i) each Ising variable is visited and a flip si — > — s; is 
proposed. The acceptance of the move is evaluated as 
in the usual Hirsch and Fye algorithm || and the new 
Green function can be computed in 0(M 2 ) operations. 

j£ j 

(ii) every independent complex variable <f> m = (<fi m ,(j) m ) 
is visited and a move <j) m — > (j> m +S^ 1 is proposed, with 
^m' 1 ^ [ — A, A] a uniformly distributed random number. 
After Fourier transforming back <j), the update requires a 
full matrix inversion (J3J), that is 0(M 3 ) operations [Io| . 
One can verify that, similarly as the original Hirsch and 
Fye algorithm, this procedure fulfills detailed balance: 
WP(W -> W) = WP(W -> W). 

Holstein phonons are bosonic fields defined by its non- 
interacting correlation function DQ(iuj m ) = —2fl /(uj' 2 n + 
fig), where fig is the phonon frequency. For the elec- 
tron conduction band that is hybridized with the impu- 
rity we shall take a semicircular density of states, thus 
A(w)/2 = t 2 cPQ (uj) = 4t 2 y/W 2 - J 1 /W 2 with Po (uj) the 
conduction electron DOS, W its half-bandwidth and t c 
is the hybridization coupling. The choice of a semicir- 
cular DOS is made because it would also be appropriate 
for our investigation of conductance through the impu- 
rity. In fact, the same A(uj) also results from attaching 
two leads to the impurity in the form of two semi-infinite 
chains of atoms with intersite hopping amplitude W/2. 
This set-up is usually used for modelling nanodevices, 
such as a molecule connected to two metallic leads. 

The DOS at the impurity p{u) = -2lm[G a (oj)] is 
obtained by a Maximum Entropy method for the ana- 
lytic continuation of the QMC calculated Green func- 
tion to the real frequency axis. We typically use upto 64 
time slices and 200,000 complete sweeps of the statistical 
fields. We consider first the simpler case of U — and 
show in Fig^systematic results for fixed A and increasing 
phonon frequency Qq. For the analysis of the results it 
is useful to recall that within the static approximation 
(i.e., taking Do(iuj m ) w Z?g(iwg)5g jTO ), which becomes 
more accurate in the small f^o limit, the integration of 
the bosons leads to an effective interaction between the 
electrons —U e ff = — 2A 2 /i7g. Since it is attractive, and 
our model is particle-hole symmetric, it will favour an 
empty or doubly occupied impurity state. Adding or re- 
moving and electron from the impurity has an energy cost 
of ~ JJ e ff/2, which is in fact the position of the two fea- 
tures in the DOS of the top panel of Fig^for the smaller 
fig = 0.333. As fig is increased and U e ff decreases, the 
DOS at low frequency has a dramatic enhancement with 
a strong quasiparticle peak emerging at u> = 0. This res- 
onance corresponds to slow charge fluctuations at the 
impurity and is the charge counterpart of the Kondo 
effect. At higher frequencies a reacher peak structure 
also evolves, where the two side peaks of the Kondo-like 
resonance correspond to the (dynamically renomalized) 
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FIG. 1: DOS for U = 0, W = 8, i c = 1. and T = 0.0357. 
Other parameters are: A = 1, flo = 0.333, 0.5, 0.75, 1., 1.5, 2., 
and A = 0.5, fio = 2 (upper to lower panel). 



Ueff/2 features, while the two additional peaks at higher 
frequency lu ~ ±(U e ff/2 + S7q) correspond to adding (or 
removing) a particle plus a phonon. When flo is further 
increased (diminishing U e ff), charge fluctuations are no 
longer prevented and the DOS moves towards that of a 
simple resonant level model with just a remanent sub- 
structure at frequencies lu ~ f^o- This basic systematic 
behavior is consistent with the T = numerical studies 
using NRG ||. 

We now turn to the evolution of the DOS as a func- 
tion of T shown in Fig|2] At high temperatures (upper 
panel) strong charge fluctuations dominate and the DOS 
is that of a resonant model with a strong and broad fea- 
ture at to = 0. In addition, two side peaks are present at 
lu ~ ±ilo due to excitations to states with an additional 
phonon. Lowering T we observe that the central feature 
splits-up into two peaks. This splitting occurs due to the 
thermal selection of the more favorable empty and doubly 
occupied states as thermal charge fluctuations decrease. 
Finally (two lower panels), bellow T\ ~ 0.0625 (anal- 
ogous to the Kondo temperature), a strong resonance 
re-emerges at to — due to the Kondo-like mechanism 
in the charge sector. ^From the Kondo model analogy, 
T\ on exp(—U e f f) therefore is rapidly suppressed with 
increasing A or decreasing f^o- 

We now switch on the interaction U and observe the 




FIG. 2: Behavior of the DOS with T for A = 1 and fi = 1. 
Other parameters are as in Fig0 
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FIG. 3: Top: DOS without Coulomb repulsion (U = 0, A = 
1). Bottom: DOS without phonons (U = 4, A = 0). Middle: 
DOS with Coulomb repulsion and phonons (U = 4, A = 1). 
The phonon frequency is Qo = 0.333, T = 0.0357 (top) and 
T = 0.0625 (middle and bottom). 



interplay between the usual Kondo effect and the vi- 
brational modes (FigEJ. Without phonons, ie, A = 
(bottom) the DOS consists of a central Kondo resonance 
and side features at ±U/2. The familiar Kondo peak is 
due to slow fluctuations between the spin-up and spin- 
down states within the charge sector of single occupation. 
Upon coupling to the vibrational modes the effective elec- 
tronic interaction of the impurity becomes dynamical 
and strongly renormalized (screened). Within the static 
approximation, its value becomes Ud yn = U — U e ff = 
U — 2 A 2 /£!o. The effective phonon interaction can even 
render the interaction attractive as in the case of su- 
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perconductivity. The effect on the impurity DOS is 
quite notable (middle), renormalizing the side peaks to 
lu ~ ±Udijn/2 and adding higher frequency peaks due 
to additional phonon excitation. A central Kondo-likc 
resonance remains, however, it is now due to a rather 
complicated many body state involving both slow spin 
and charge fluctuations. 




FIG. 4: Conductance (a. u.) at T — 0.0625 and A = 1 as 
a function of the gate potential Vg. Left: U — and f2o = 
0.333, 1,2 (bottom to top). Right: U = 4 and fi = 0.333. 
For comparison, in thin line we show p(u>) for V g = 0. 

We may now turn to the discussion of the consequences 
of the above mentioned effects on the conducting proper- 
ties of this quantum impurity model that should be rele- 
vant for nanoscopic devices such as a molecule attached 
to leads with both electronic and vibrational degrees of 
freedom. 

One of the usual experimental setups to study trans- 
port phenomena in nanostructures U 0, 0] i consists in 
applying a bias potential between the left and the right 
sides of the impurity (molecule) , and then measuring the 
current through it. In addition a gate voltage V g might 
be applied at the impurity site, and the conductance is 
obtained as a function of V g . 

In the limit of very small bias between the left and 
right leads the conductance can be calculated from Q 




where f(u) = 1/(1 + e^). Results at low T are shown 
in Fig^J In this parameter regime, the behavior of G is 
essentially given by the value of p(0). For U = (left) 
a maximum in the conductance is observed at V g = 0, 
and its value is strongly enhanced as U e ff is decreased. 
The maximal value of the conductance is controlled by 
the spectral strength of the Kondo-like resonance that is 
rapidly suppressed by U e f / • While the V g = case is un- 
derstood in terms of the analogy to the usual Kondo case, 
we find that upon application of a gate potential there is 
a very different behavior with respect to the case of con- 
ductance through a magnetic impurity. In the latter case, 



at low T, the conductance displays peaks at ±t//2 and 
remains sizable in between |4J . This is because the Kondo 
peak survives the departure form the particle-hole sym- 
metric situation, since V g does not break the spin doublet 
at the impurity and the slow spin fluctuation mechanism 
remains operational. In stark contrast, the present case 
shows a rather narrow peak of the conductance centered 
at V g — and this is due to the fact that a finite V g imme- 
diately destroys the empty - doubly occupied doublet at 
the impurity, and the validity of the Kondo-like analogy. 
At low temperatures, the probability of the participa- 
tion of those states differs by a factor ~ exp(—V g /T). 
Switching on U (right), we find a new surprise, namely 
that even for Ud yn < (ie, the phonon mediated in- 
teraction dominates) the conductance remains high in a 
neighborhood of V g = and develops a two peak struc- 
ture at V g « This behavior is consistent with 
a Kondo-likc character that remains present even away 
from particle hole symmetry. The systematic investiga- 
tion of the conductance in the whole parameter regime 
is left for future investigations. 

In conclusion, we have introduced a new QMC algo- 
rithm for the numerical solution of a general class of 
quantum impurities models with fermionic and bosonic 
degrees of freedom. We considered the specific case of a 
SIAM with a Holstein phonon that served both to bench- 
mark the quality of our method and to demonstrate new 
effects in the impurity conductance. This model and fur- 
ther generalizations might be of great help for the in- 
vestigation of conduction properties through molecular 
nanodevices. 
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